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Abstract 

Bouncing jets are fascinating phenomenons occurring under certain conditions when a jet impinges on a free surface. 
This effect is observed when the fluid is Newtonian and the jet falls in a bath undergoing a solid motion. It occurs also 
for non-Newtonian fluids when the jets falls in a vessel at rest containing the same fluid. We investigate numerically the 
impact of the experimental setting and the rheological properties of the fluid on the onset of the bouncing phenomenon. 
Our investigations show that the occurrence of a thin lubricating layer of air separating the jet and the rest of the liquid 
is a key factor for the bouncing of the jet to happen. The numerical technique that is used consists of a projection 
method for the Navier-Stokes system coupled with a level set formulation for the representation of the interface. The 
space approximation is done with adaptive finite elements. Adaptive refinement is shown to be very important to 
capture the thin layer of air that is responsible for the bouncing. 

Keywords: Bouncing Jet, Kaye effect, Entropy Viscosity, Level set, Projection Method, Shear-thinning viscosity, 
Adaptive Finite Elements 


1. Introduction 

The ability of a jet of fluid to bounce on a free surface has been observed in different contexts. Thrasher et al. Oil 
designed an experiment where a jet of Newtonian fluid falls into a rotating vessel filled with the same fluid. The authors 
investigated the conditions under which the jet bounces: nature of the fluid, jet diameter, jet and bath velocities. We 
refer to Lockhart l26l for experimental movies illustrating this effect. Bouncing can also be observed in a stationary 



Figure 1: Experimental Observation of the Kaye Effect. (Left) the fluid start to buckle producing a heap; (middle) A 
stream of liquid suddenly leaps outside the heap; (right) Fully developed Kaye effect. 

vessel provided the jet is composed of a non-newtonian fluid. During the pouring process, a small heap of fluid forms 
and the jet occasionally leaps upward from the heap; see Figure [T] This is the so-called Kaye effect as first observed 
by Kaye El in 1963. About 13 years after this phenomenon was first mentioned in the literature, Collyer and Fischer 
m revisited the experiment and suggested that the ability for the fluid to exhibit shear-thinning viscosity and elastic 
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behavior are key ingredients for the bouncing to occur. Additional laboratory experiments performed by Versluis et al. 
(Ml and Binder and Landig 0 lead the authors of each team to propose a list of properties that the fluid should have 
for the Kaye effect to occur. The conclusions of these two papers disagree on the requirement that the fluid be elastic 
and on the nature of the thin layer that separates the heap and the outgoing jet. It is argued in Versluis et al. f34l that 
elastic properties are not necessary and that the thin layer is a shear layer, whereas it is argued in Binder and Landig 
0 that the fluid should have elastic properties and that the thin layer separating the heap from the bouncing jet is a 
layer of air. Recent experiments reported in Lee et al. f25l using a high speed camera unambiguously show that the 
jet slides on a lubricating layer of air. For completeness, we also refer to Ochoa et al. 1771 for a thorough discussion 
on the “stable” Kaye effect, where the jet falls against an inclined surface. 

The objective of the present paper is to numerically revisit the Kaye effect. Our key finding is that a thin air 
layer is always present between the bouncing jet and the rest of the fluid whether the fluid is Newtonian or not. Our 
numerical experiments suggest that the critical parameter for bouncing to occur is that the properties of fluid and the 
experimental conditions be such that a stable layer of air separating the jet and the ambiant fluid can appear. This 
condition is met by setting the bath in motion for Newtonian fluids; it can also be met if the bath is stationary provided 
the fluid is non-Newtonian and has shear-thinning viscosity. The numerical code that we use is based on a modeling 
the fluid flows by the incompressible Navier-Stokes equations supplemented with a surface tension mechanism. The 
shear-thinning viscosity of the non-Newtonian fluids is assumed to follow a model by Cross E). Elastic behaviors are 
not modeled. The numerical approximation of the resulting two-phase flow model is based on two solvers: one solving 
the Navier-Stokes equations assuming that the fluid/air distribution is given, the other keeping track of the motion of 
the interface assuming the transport velocity is given. The Navier-Stokes solver is based on projection method by 
Chorin El and Temam l30l using a second-order backward differentiation formula for the time discretization and 
finite elements for the space approximation. The transport solver is based on a level-set technique in the spirit of Osher 
and Sethian l28ll to represent the liquid/air interface. The level set is approximated in space by using finite elements 
and the time stepping is done by using a third order explicit Runge-Kutta technique. 

The paper is organized as follows. The mathematical model is presented in Section [2] The numerical techniques 
to solve the Navier-Stokes equations and the transport equation for the level set function are described in Section [3] 
Various validation tests of the numerical algorithms and comparisons with classical benchmark problems are reported 
in Section]?] Finally, we report numerical evidences of Newtonian and non-Newtonian bouncing jets in Section]?] 

2. The Mathematical Model 

This section presents the mathematical models adopted to describe non-mixing two phase fluid flows with capillary 
forces. Each fluid is assumed to be incompressible. 

2.7. Two Phase Flow System 

Let A c R d (d = 2,3) be an open and bounded computational domain with Lipschitz boundary DA and let 
[0, T ] be the computational time interval, T > 0. The cavity A is filled with two non-mixing fluids undergoing some 
time-dependent motion, say fluid 1 and fluid 2. We denote by D + and Q~ the open subsets of the space-time domain 
Ax[0, T] occupied by fluid 1 and fluid 2, respectively. We denote by p + , p + , p“, P _ the density and dynamical viscosity 
of each fluid, respectively. The interface between the two fluids in the space-time domain is denoted 2 := <9Q + Pi d£2~, 
and the normal to 2, oriented from Q + to Q", is denoted n^. The two space-time components of the vector field 
nx : 2 —■> are denoted n and n T , respectively. It is also useful to define Q ± (0 := f 1 ± n (Ax{r}); i.e., the sets 

Q + (0 and Q~(t) are the regions occupied by fluid 1 and fluid 2 at time t, respectively. We also introduce the interface 
Z(t) = dQ + (t) P dfl~(t ); note that n, as defined above, is the unit normal of 2(t) and it is oriented from £2 + (7) to 
£l~(t). To facilitate the modeling, we define global density and dynamical viscosity functions p,p : Ax[0, T] —> R by 
setting p(x, t) = p ± and p(x, t) = p ± if (x, t) e Tl ± . The fluid velocity field u : Ax[0, T] —> R d , henceforth assumed 
to be continuous across 2, and the pressure p : Ax[0, T] —> R are defined globally and solve the incompressible 
Navier-Stokes equations in the distribution sense in the space-time domain 

p|-^u + ii-Vuj - 2div(p + Vp - S^cnai = pg in Ax(0, T], (la) 

div(u) = 0 in Ax(0, T], (lb) 
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where V^u := |(Vu + Vu r ) is the strain rate tensor, g is the gravity field, and S^cr/ai is a singular measure modeling 
the surface tension acting on X(7). The distribution 6% is the Dirac measure supported on S, the function k : S —> M is 
the total curvature of S(0 (sum of the principal curvatures) and <x is the surface tension coefficient. 

The system Q is supplemented with initial and boundary conditions. The initial condition is u(x, 0) = uo(x) for 
all x e A, where Uo is assumed to be a smooth divergence-free velocity field. The boundary dA is decomposed into 
three non-overlapping components <9 A := T D U T N U Y N with T D nT N = 0, T^nT^ = 0, Ts n T# = 0. Time- 
dependent decomposition could be considered but are not described here to avoid unnecessary technicalities. Given 
f/v : T n —> R d and f D : T D —> R d , we require that 
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Co 
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on T n x( 0, T], 

(2a) 

U = f D , 

on Tz)X(0, T], 

(2b) 

uv = 0, ((2^V s u - pl'j v) xv = 0, 

on T s x(0, T], 

(2c) 


where v is the outward pointing unit normal on dA and I is the dxd identity matrix. For simplicity, we assume that 
the (d - l)-measures of T N U Ts and T# are each strictly positive; otherwise extra constraints either on the velocity or 
on the pressure must be enforced. 

Note that we could have formulated the conservation of momentum without invoking the singular measure mod¬ 
eling the surface tension by saying that •DU 1 holds in Q + and Q (without the singular measure) and by additionally 
requiring that 

[u] = 0 and - /?/] n = cr/cn on S(0, W 6 (0, T], (3) 

where [.] denotes the jump across S(0 defined by [v](x) := linto- 9y _> x v(y)-limQ+ 9z _> x v(z), i.e., [v](x) = v(x“)-v(x + ), 
for all x e S(0 and all v : A —> R d or v : A —> R dxd . 

The interface Z(t) is assumed to be transported by the fluid particles. More precisely, let <9“(Ax[0, T]) be the part 
of the boundary of the space-time domain where the characteristics generated by the field (u, 1) enter, i.e., 

<9“(Ax[0, T]) := Ax{0} U {(x, t) e <9Ax[0, T] | u(x, t)-v < 0}. (4) 

We then define {x(P, t) e A, t e [ t P , T], (P, tp) e <9“(Ax[0, T])} to be the family of the characteristics generated by 
the velocity field u, i.e., f x(P, t) = u(x(P, t), t) with x(P, tp) = P, (P, tp) e <9“(Ax[0, T]) where tp is the time when the 
characteristics enters the space-time domain Ax[0, T] at P. Let us now denote by So := 2 n <9“(Ax[0, T ]) the location 
of the interface at the inflow boundary of the space-time domain, then we are going to assume in the entire paper that 
the velocity field is smooth enough so that the following property holds 

Vx e S(£), 3!(P, t P ) e S 0 , x = x(P, t), t > t P . (5) 


2.2. Eulerian Representation of the Free Boundary Interface 

A level set technique is used to keep track of the position of the time-dependent interface S(£), see for instance 
Osher and Sethian [28]. This method is recalled in Section [2.2.1 A typical problem arising when using a level-set 
method to describe interfaces is to guarantee the non-degeneracy of the representation. We discuss a reinitialization 


technique overcoming this issue in Section 2.2.2 


2.2.7. Level-Set Representation 

Let us define the so-called level function 0(x, t) : Ax[0, T] —> R so that 

-f-0 + u-V0 = 0 in Ax(0, T], 0(P, tp) = 0o(P, tp) on <9“(Ax[0, T]), 

dt 

where we assume that 0o is a smooth function satisfying the following properties: 


dQf n d~ (Ax[0, T]) = {(P, tp) e 8~ (Ax[0, T]) \ ± 0 O (P, t P ) > 0}. 
So := S n <T (Ax[0, T]) = {(P, tp) e d~ (Ax[0, T]) | 0 O (P, fp) = 0}. 


( 6 ) 
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(7) 

( 8 ) 






Note that this definition implies that 2 0 is the 0-level set of 0 q. Upon introducing if/(F,t) := 0 (x(P, t), t) for all 
(P, fc) E d~(Ax[0, T]), t > tp, the definition of 0 together with the definition of the characteristics implies that 
d t if/(F, t) = 0; thereby proving that 0 (x(P, t), t) = 0 (P, t) = i//( P, t v ) = 0 (x(P, tp), tp) = 0 (P, tp) = 0o(P, tp). This means 
that the value of 0 (x(P, t), t) along the trajectory {x(P, tp ), t e [tp , T]} is constant; in particular the sign of 0 (x(P, t), t) 
does not change. The leads us to adopt the following alternative definitions for D* and 2(0: 

Or := {(x, t) e Ax[0, T] \ ± 0 (x, t) > 0} 

2(0 := {x e A | 0 (x, 0 = 0}. 

The above characterizations will be used in the rest of the paper. 

To simplify the presentation we are going to assume in the rest of the paper that there is T L c dA such that 
<9“(Ax(0, T)) = T L x(0, T), i.e., the inflow boundary for the level-set equation is time-independent. We then set 
0 L (P, 0 = 0o(P, 0 for all F ET L ,t E [0, T], and we abuse the notation by setting 0 O (P) = 0o(P, 0) for all p € A. 

2.2.2. Reinitialization and cut-off function 

A typical issue when dealing with level-set representation of interfaces is to guarantee that the manifold {x e 
A | 0 (x, 0 = 0} is (<d - l)-dimensional, i.e., we want to make sure that ||V 0||^2 > 0 in every small neighborhood of the 
zero level-set, where || • \\p is the Euclidean norm in W*. In order to achieve this objective, we implement an “on the 
fly” reinitialization algorithm proposed in Ville et al. E3, which consists of replacing ([ 6 ]) by 

= Tsign(0^) - ||V0^||^), (10) 

where A,(3 > 0 are parameters yet to be defined, and ^(-), sign(-) are defined by 

2 ( 1, z>0 

, sign(z) := < -1, z<0 (11) 

{ 0, \z\ = 0. 

The rational for the new definition ( fTT)] ) is that the presence of the sign function in the right-hand side implies that 
the 0 -level set of is the same as that of 0 ; i.e., the characterizations of Q* and X(t) are unchanged, see 
Moreover, assuming that u is locally the velocity of a solid motion and upon setting ^(P, t) = (p^(x(F, t )), we have 
j t if/ = Tsign( 0 r )((l - (|) 2 ) - |V 0 j). Assuming that this eikonal equation has a steady state solution and denoting by 
2 oo the 0 -level set of this steady state solution, the behavior of in the vicinity of 2 ^ is ^(p) ~ (3 tanh ^ dlst ^°°’ p ) ), 
since the solution to the following ODE : y'(z) = 1 - (y(z)//3 ) 2 is 

T(z):=p tanh . (12) 

Note that if/ is close to if/oo if A is very large. In conclusion the solution of (T 0 | ) is such that 0 ^(x, t) « (3 tanh ( dlst( ^ (r),x) ) 
if A is large enough, i.e., ||V 0 ||^> ~ 1 in any small neighborhood of the zero level-set. We are going to abuse the 
notation in the rest of the paper by dropping the indices A, (3 and by using 0 instead of 0 ^. 




(9a) 

(9b) 


3. Numerical Method 


We disc uss the approximations of the level-set equation in Section 3.1 and that of the Navier-Stokes system in 


Section 


3.2 


We consider a mesh family, [Th)h> o» and we assume that each mesh Th is a subdivision of A made of 
disjoint elements K , i.e., rectangles when d = 2 or cuboids when d = 3. We denote by S(Th) the collection of interfaces 
and boundary faces (edge when d - 2 and faces when d - 3). Each subdivision is assumed to exactly approximate 
the computational domain, i.e., A = U K ^ h K, and to be consistent with the decomposition of the boundary, i.e., there 
exists Si(T h ) c S(Th) for I e {D,N,S,L} such that Tj = U fes^F- The diameter of an element K e Th is denoted 
by hx\ hx := max X; y GjS : |x - y|. The mesh family {Th)h >o is assumed to be shape regular in the sense of Ciarlet. For 
any integer k > 1 and any ^ e 7^, we denote by Q^(^T) the space of scalar-valued multivariate polynomials over K of 
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partial degree at most k. The vector-valued version of Q k (K) is denoted Q k (K). The index h is dropped in the rest of 
the paper and we write T instead of Th when the context is unambiguous. 

Regarding the time discretization, given an integer N > 2, we define a partition of the time interval 0 =: t° <t l < 
... <t N :=T and denote 5f :=t n - t n ~ l and t n+ i := t n + . 


3.1. Numerical Approximation of the Level-Set System 

The continuous finite element method used for the space approximation of the level-set equation © is described 


in Section 3.1.1 We present in Sections 3.1.2| an entropy-residual technique that has the advantage of avoiding 


the spurious oscillations that would otherwise be generated by using an un-stabilized Galerkin technique. The time 
stepping is done by using an explicit strong stability preserving (SSP) Runge-Kutta 3 (RK3) scheme as explained in 
Section [3X3] 


3.1.1. Approximation in Space 

The space approximation of ®(-, t) of the level set function 0(-, t) solution to ( fTO] ) is done by using continuous, 
piecewise linear polynomials subordinate to the subdivision T. The associated finite element spaces are defined by 

W (T) := {W e C°(A; R) | W\ K e Q\K), VK e T), (13) 

WoCT) := {W g C°(A; R) | W = 0 on r L , W\ K e Q l (K), VK e T), (14) 

W l (T) := {W g C°(A; R) | W = <D L on T L , W\ K e Q\K), g T), (15) 

where O l (-, t) is a piecewise linear approximation of the inflow data 0 L (-, t). Assuming that the velocity field u : 
Ax[0, T ] — > R d is known, the Galerkin approximation of ( fTO] ) is formulated as follows: Given and 0(-, 0) = O 0 , 
where ®o £ W(7~) is an approximation of the initial condition 0o, find ® e CHtO, T]; Wl( 7~)) such that 

f 2,d w = _ f u-V<DW+ f Asign(<D)(^(<D)-||V<D|^)W, VW e W 0 (7"). (16) 

Ja at J A J A 

It is well known that the solution to the above system exhibit spurious oscillations in the regions where ||VO ||^2 is 
large. We address this issue in the next section. 


3.1.2. Entropy Residual Stabilization 

We describe in this section an entropy viscosity technique to stabilize the Galerkin formulation ©• This method 
has been introduced in Guermond et al. ED and we refer to Bonito et al. m for a mathematical discussion on its sta¬ 
bility properties. To motivate the discussion, we refer to the panel (b) of Figure [2] showing the Galerkin approximation 
of the characteristic function of the unit disk, initially centered at (|, 0), after one rotation about the origin. 

The spurious oscillations are avoided by augmenting ( fl6| ) with an artificial viscosity term where the viscosity is 
localized and chosen to be proportional to an entropy residual. To describe the method and define an appropriate local 
“viscosity”, we recall that the following holds in the distribution sense for any E e C X (A; M): 

+ u-VE(0) - 4sign(0)(£(<D) - ||V4>\\ ( i) E'($) = 0, 
ot 

Consequently, it is reasonable to expect that the semi-discrete entropy residual 

R™(A\ u) := --Em + U(r)-V£(0) - isign(O) (Qm - ||V®||W E'(0)b 
Ot 

is a reliable indicator of the regularity of f. This quantity should be of the order of the consistency error in the regions 
where 0 is smooth and it should be large in the region where the PDE is not well solved. In our computations, we 
have chosen 

Etf) = W, (17) 
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(C) (d) 


Figure 2: Graph of the level-set function in the computational domain A = (-1, l) 2 using the velocity given by e 
05): (a) Initial data: 0o is the characteristic function of the disk of radius 1 centered at (^, 0); (b) No stabilization; (c) 
First-order stabilization with Cun = 0.1; (d) Entropy viscosity stabilization with Cun = 0.1 and Ce nt = 0.1. Observe 
the spurious oscillations in panel (b) when no stabilization is applied. Both, the first-order and the entropy viscosity 
solution are free of oscillations, the latter being clearly more accurate. 


The local so-called entropy viscosity is defined for any K e T by 

||/? Ent (<f>,u)|b 


A4 nt (®,u) := C EDt h 2 K - 


(K) 




( 18 ) 


where CEnt is an absolute constant. In the regions where 0 is discontinuous (or has a very sharp gradient), the entropy 
viscosity as defined above may be too large and thereby introduce too much diffusion, which in turn may severely 
limit the CFL number when using an explicit time stepping. In this case a linear first-order viscosity is turned on 
instead 


/4' n (®, u) = C Lin fe(u + dsigm®)- 


V® 




(19) 


IIVOII^ 

where Cun is an absolute constant. The justification for the definition of the local speed that is used to define the 
viscosity in © is that ([TO]) can be re-written ^0 + w-V0 = Tsign(0)^(0), where w = u + Tsign(0) ||V ^ 2 . Combining 

the two viscosities yield the artificial viscosity yu stab : Ax[0, T] —> M defined on each K e T by 

/i Stab (<F,u)|* := min0uL in (d),u),4 nt ((I),u)). 

Going back to the space discretization, we modify (\6^ as follows: Look for ® 6 C^tO, T]; Wl(7~)) so that 

' Stab (0,u)V0VlT 


f Uw = - f u-VO W + f Tsign(O) (^(O) - ||VO||^) W - f 
J a °t J A J A J A 

for all W e Wo (T) and d>(-, 0) = O 0 . 

3.1.3. Approximation in Time 

Before introducing the time discretization we re-write pT| ) as follows: 

f |-®w= f L(®,u)W- f jU Stab (®,u)V® VW, 

Ja at J A J A 




( 20 ) 


( 21 ) 


VVF e Wq(T, t). 


( 22 ) 


where L(0, u) := -u + Tsign(O) (^(O) - ||VO||^). Then we approximate time in the above nonlinear system of ODES 
by using an explicit RK3 strong stability preserving (SSP) scheme, e.g. see Gottlieb et al. (HI, Shu and Osher 


6 







for more details on SSP methods. We denote by the approximation of ®(-, t k ), 0 < k < N. Then, the time stepping 
proceeds as follows: Given ® w , compute O 1 , ® 2 ,0 3 e W(T) and O w+1 e Wl(7~) so that 


(0 

(«) 

(m) 

and 


J <D (1) W = J (<D" + (5r" +1 L(®",u' I ))w, VffeWfT) 

J 0 (2) W = J > |^<D" + ^<D (1) + ^r' !+1 L(<D (1) ,u" +1 )jw-J" jU Stab -" +1 VO (1) -Viy, 
J ® (3) W = J |i<D" + ^<D (2) + ^f)Y ,+l L(O l2l ,u' ,+ 5)j W - Ju Stab ’ n+ Wo (2) -VW, 


<D n+1 (x,t) = 


<D (3 ) (x) 


when x£T L , 


<D L (x, t n+L ) when x eT l . 


VW 6 W (T) (23) 
VW e W(T) 


(24) 


Using ( [20] ), the viscosities ^ Stab ’ n+1 and ^ Stab ’ w+ 2 are defined to be equal to yu Stab (® (1) , u n+1 ) and yU Stab ( 0 (2) ,u w+ 2 ), 
respectively, where the residual R Ent in the definition ( p~8j ) of p Ent is evaluated as follows: 


fl(® (1) , u" +1 ) * £((I)I If* ) + (u(/" +1 )-V<D (1) - Asign(<D (1) ) (@(<D (1) ) - ||V<D (1) ||^)) £'(<D (1) ) (25) 

/?(0 (2) ,u' t+ 5) * E(<1)< i^. ~ + f (fl> ^ + (u(?" +2 )-V1> (2) - isign(® (2) ) (@(® (2) ) - ||V<D (2) || f 2 )) £'(® (2) ). (26) 


Remark 1 (“On the Fly” Stabilization). Notice that following Bonito et al. Ml/, no viscosity is added in the computa¬ 
tion of& l \ In particular, the viscosities used within the time interval [ t n , t n+1 ] only depend on the values o/O and u 
on the same time interval. 


Remark 2 (Stability and Convergence). We expect the scheme to be stable under the following Courant-Friedrichs- 
Lewy (CFL) condition: 

n +1 - ^ ...... h K 


5f + 


< CcFL m i n 

KeT 


\u(P 


n+l\ 


(27) 


) + 2sign(O w ) ||von|| ^ ||l-(^) 

for some sufficiently small but positive constant Ccfl independent of A, T, St, ® and u. Refer for instance to Bonito 
et al. Ml/ for further details on the CFL condition. Moreover, it seems reasonable to expect that only the entropy 
viscosity is active in the regions where f is smooth; as a result, the CFL condition implies that the first-order ap¬ 
proximation of the time derivative in the evaluation of the entropy residuals, ([25])-([26]), does not affect the overall 
third-order approximation thanks to the h 2 K factor present in the definition of the viscosity CD- This conjecture is 


confirmed numerically in Section 4.2 


3.2. Numerical Approximation of The Navier-Stokes System 

The space approximation of the velocity and pressure in the Navier-Stokes equations is done by using Taylor-Hood 
finite elements. The time discretization is done by using the second-order backward differentiation formula (BDF2). 
An incremental rotational pressure correction scheme is adopted to uncouple the velocity and the pressure. We refer 
to Guermond and Shen (20 ] for the convergence analysis of the metod and to Guermond et al. l2TTl for a review on 
projection methods. 


3.2.1. The Space Discretization 

Let F d be a continuous, piecewise quadratic approximation of f# on T#. The finite element discretization of the 
velocity and the pressure is done by using the following linear and affine spaces: 

V 0 (T) := {V € C°(A; R d ) | V|r D = 0, V-v|r s = 0, V|* € Q 2 (Z0, VK e T), (28) 

VdCD := {V € C°(A; R d ) | V|r D = F 0 , V-v|r s = 0, \\ K € Q 2 (K), VK e T}, (29) 

M(T) := [Q 6 C°(A) —> R | Q\ K € Q 1 (K), VK e T), (30) 
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Upon setting U(0) = Uo, where Uo is a continuous, piecewise quadratic approximation of the initial velocity Uo in 
V£>(7~), the semi-discrete formulation of ([]} consists of looking for U e C^tO, T];V D (T)) and P e C°([0, r];M(7~) 
such that the following holds for every t e (0, T ]: 


£ p|^U + (U-V)u| V + 2 £ n (V S U:V S 

= f pg-V+ f f N -V+ f (tk n-V, 

J a Jr N Jm) 


V)-J* PdivV + £ QdivU 
V(V, Q) e VoCDxMCT), 


(31) 


3.2.2. Time discretization 

The time discretization of ( [3T] ) is done by using the second-order backward differentiation formula (BDF2) and an 
incremental pressure correction scheme in rotational form introduced and studied in Guermond and Salgado m to 
decouple the velocity and the pressure. In addition to the approximation of the initial condition U° mentioned above, 
the algorithm requires an approximation P 0 e M(T) of the initial pressure p( 0). We denote by U n , P n and x ¥ n the 
approximations of U(-, t n ) and P(-,t n ) and the pressure correction, respectively. 

The initialization step consists of setting: U° = Uo, P~ l = P° = Po and T* 0 = 0. Then, given a new time step St n+l , 
given U n e V D (T), ¥ n e M(7“) and P n e M(7“), and assuming for the time being that p(t n+1 ), p(t n+1 ) and £(7 n+1 ) are 
also known (see ^3.4[ ), the new fields XJ n+l e V^(7“), X F W+1 e M(7“) and P n+l e M(7“) are computed in three steps: 

\±\ Velocity Prediction: Find U n+1 e Vd( 7“) such that 


f p(^ +1 )-^-U" +1 -V + 2 f p(^ +1 )(v s (U" +1 ):V s v) + 5 r (U" +1 ,V) 

J A at Ja 


= - f p(f" +1 )((U")*-VU")-V + f(P n + - I'P" _1 )div(V) 

J A J A 3 3 


■ f p n+l g n+1 + f ftf +1 -V + f cr/( n+1 n" +1 -V, VY € Vo(7~), 

J a Jr N Jx(t n+1 ) 


(32) 


where 

• the BDF2 approximation of the time derivative with variable time stepping is given by 


4lJ(x,f" +1 ) - 

dt St n+l 


BDF2 TT^+l 


U (x) := 


1 (\+ 2//„ + 1 


5t n+l ^ 1 + Tj n+ i 


U" +1 (x) - (1 + ? 7 „ + i)U"(x) + 


1 + Vn+l 


U" _1 (x)|, 


m +1 


St‘ 


-/+1 


St 1 


• the extrapolated velocity (U n )* is defined by (U n )* := \3 n + if(U n - X3 n 1 ); 

• as discussed in Bonito et al. 0, the bilinear form Sr • V(7~)xV(7~) —> M is added to control the divergence of 
the velocity and to cope with variable time stepping and open boundary conditions 


S r (W,V):=C stab y] f(p" +1 +p' ,+1 |feW|| L » m )(V-W)(V-V), 

KeT ^ K 


where C sta b is an absolute constant. 

PH Pressure Correction Step: The pressure increment x F n+1 e M(7~) is determined by solving 
£ VT" +1 V0 = - 3mm 2 ^ff" +l > £ div(U n+1 )<2, V2 € M(T). 

0 Pressure Update: The pressure P n+l e M(7~) is obtained by solving 

£ p n+l Q = £ [p n + >P" +1 )Q-minyu(f" +1 ) £ div(U" +1 ) Q, VQ e M(T). 


(33) 


(34) 


(35) 
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For stability purposes, we restrict the space and time discretization parameters to satisfy a CFL condition 


<S/" +1 


< C’cFI. 


1 mm KeTh h K 

2 ||U»|| L ~ (A) ’ 


( 36 ) 


where Cqfl is the same constant appearing in ( |27] l and j m\n KeTh Iik is the minimum distance between two Lagrange 
nodes using Q 2 elements. 


3.3. The Surface Tension 

This section describes the approximation of the curvature term JT crK n+l (n w+1 -V) appearing in the first step of 
the projection method ( [32] ). We follow the method proposed by Bansch [4] (see also Hysing et al. [23]) using the work 
of Dziuk and Elliott fl4l . This approach is based on the following representation of total curvature: 


Kn = Vx-V^Id^, 


where Id^ is the identity mapping on X and, given any extension v of v in a neighborhood of X, the tangential gradient 
of v : X —> is defined by V s v := Vv| s (I - n ® n), see e.g. Gilbarg and Trudinger 1171 . Multipliying the above 
identity by a test function V and integrating by parts over X yields 


| cr/cn V = - | crV^Ids : V S V + | crd T Id 2 -V, 

Js Js J<9£ 


(37) 


where r is the co-normal to X and d T is the derivative in the co-normal direction. In the present context the integral 
on <9X vanishes since either X is a closed manifold or d T Id 2 = 0. This identity together with the first-order prediction 
Id S ( ? n+i) « Ids(^) + St n+l \J n+l of the interface evolution 0 gives a semi-implicit representation of the total curvature 



o-K n+l (n w+1 -V) = - 



V jxft+i ) Id^+i): V 2(*n+i) V 
V z(t n+] )Ids(^+i): V) V - St n+l 



^V^ +1) U W+1 :V^ +1) V. 


One key benefit of this representation comes from the additional stabilizing term j£ j crV £ (^+i)U n+1 : V^+qV, which 
we keep implicit in ( [32] ). The technicalities regarding the approximation of this integral using the level set representa¬ 
tion are detailed in section [3A2] see ( [45] ). 


3.4. The Coupled System 

The two solvers described in Sections [3T] and [T2] above are sequentially combined. The flowchart of the resulting 
free boundary flow solver is shown in Figure [3] The remaining subsections of jj |3.4| detail the coupling between the 
two solvers and other implementation technicalities. 
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3.4.1. Data for the Lev el-Set Solver <[33)-{25 

Following ( [23]) , given <D W , an approximation of f(t n ), the approximation <D" +1 of 0(^ +1 ) requires the velocities 
u(P), u(P + 2 ) and u (t n+1 ). To avoid an implicit coupling between the level-set solver and the Navier-Stokes solver, 
these quantities are replaced by second-order extrapolations using U" and ll n ~ l : 

i 1 St n+l St n+l 

u(/")~U\ u(/' ,+ i)»U" + --—(If-lT -1 ), and U" +1 « U" + —— (If - U” -1 ). 

2 ot n ot n 

The sign(.) function that is used in the right-hand side of (TO] to make sure that ||V0||^> is close to 1 in a small 
neighborhood of £ (i.e., on the fly reinitialization, see discussion following fTT] ) is redefined and replaced by: 

( +1 if s >/3tanh(Cs), 

sign h (s) = l -1 if s < -y?tanh(C 5 ), (38) 

{ 0 otherwise, 

where Cs is an absolute constant. The thresholding in the above definition of the approximate sign function gives 
sign /2 (0(x, t)) = ±1 whenever > tanh(Q), which is compatible with the behavior « tanh that is 

expected for the level-set function, see (T2] . 

The parameter A (“reinitialization relative speed” in the language of Ville et al. (35]) is defined for t e (t n , t n+l ] by 

3 = C x I|U1l~ ( a)> KeT, (39) 

where Cx is an absolute constant. This definition is motivated by the CFL condition ([27]). We refer to Ville et al. (35] 
for further details. 


3.4.2. Data for the Navier-Stokes Solver 

The definitions of the fields \J n+l and P n+l in ([32])—([35]) invoke the values of the density field p n+l and viscosity 
field ji n+l . Once f n+l is computed, these quantities are evaluated by using the following definitions 


n+ i + 1 + W>" +1 ) _1 

P =P -~-+ P - 


,n +1 


+ 1 + Hh((^ n+l ) _1 

: P -o-+ P “ 


2 


(40) 

(41) 


2 2 

where, p ± ,p ± are the density/viscosity in Q*, and Hh(.) an approximation of the Heaviside function defined as follows: 


H h (s) 


1, 

- 1 , 


if s > /?tanh(C#), 
if s < -/Hanh(Ctf), 

otherwise, 


(42) 


ygtanh(C^)’ 

where C H is an absolute constant, as suggested in Ville et al. (35] . Similarly to what we have done to approximate the 
sign function, the above regularization is compatible with the behavior « tanh ( dlst( Vd,x) ^ j s expected for the 
level-set function. 

The approximation of the surface tension term in (32] is done by following Hysing et al. [23]. Let e > 0, and 
consider the piecewise linear regularized Dirac measure supported on £(P +1 ), say S e , defined by 


S £ (x) : 


( I ^ _ dist(£,x) ^ 


if |dist(2,x)| < 6, 
otherwise. 


(43) 


To account for the fact that we do n ot have access to the distance to the interface 2 but rather to an approximation 


0 f^tanh(^f^ 
instead: 


), see Section 


2.2.2 


6 e (x, O) 


we rescale S e (x) as suggested in Engquist et al. m, Tornberg (331 and consider 

(44) 


= | l(l- £ F)llVO(x)|| f2) l®(x)|<e:=egf 


6 

o, 


otherwise, 
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where || • W& is the ^-norm in R d . In practice, we chose e = /3tanh(C#) to be consistent with the approximation of 
the Heaviside function, see ( [42] ). Using this approximate Dirac measure, the approximation of the surface tension 
discussed in Section 1331 becomes 


f cTK n+l (n w+1 -V) « - 

JYj(t n+l ) 


f 

I 


V xu'i - 1 )Idv ( ,»+1 } : Vv (( „+1 j V - St" 
(V 0 „ + iId A :V fl ^iV)tf e (.,<D" +1 ) - 5f 


'X 


*'X( 


<TV S( ,„ +1) U n+1 :V E(( „ +1) V 
crV^n+i U" +1 : V) S £ {., <D" +1 ), 


(45) 


where ld\ is the identity mapping on A and 


V<D«+iIdA — I ~ 


VO w+1 0 VO n+1 
l|VO w+1 ||^ 2 


V 0 „ + iV := VV 


I- 


VO w+1 0 vo w+1 " 
l|VO" + 1 |l * 2 , 


VV € Vd(T). 


Note that the above definitions correspond to approximating the normal vector n on l<(t n+l ) by - 2 


3.5. Adaptive Mesh Refinement 

The experiments reported in Thrasher et al. ED, Binder and Landig s, Lee et al. |[25l suggest that a critical 
feature of bouncing jets is the occurrence of a thin layer between the jet and the rest of the fluid. These observations 
have led us to adopt a mesh refinement technique to describe accurately this thin layer. The refinement strategy is 
designed to increase the mesh resolution around the zero-level set of <D. More precisely, a cell K is refined if its 
generation count (the number of times a cell from the initial subdivision has been refined to produce the current cell) 
is smaller than a given number R max and if 


mx K ,t)\</3tonh(C R ), (46) 

where x K is the barycenter of K and C R is an absolute constant. The purpose of the parameter R max is to control the 
total number of cells. The parameter Cr controls the distance to the interface below which refinement occurs. Note 
that ( [46] ) is compatible with ( [38] ) and ( [42]) . The subdivisions are done with at most one hanging node per face, see 
e.g. Bonito and Nochetto Section 6.3]. A cell is coarsen if it satisfies the following three conditions: its generation 
count is positive; 

|O(x*,0l>/Hanh(C c ), (47) 

where Cc > C R is an absolute constant; and if once coarsened, the resulting subdivision does not have more than one 
hanging node per face. We refer to the documentation of the deal.II library for further details, Bangerth et al. 0. 

4. Numerical Validations 

The algorithm presented in the previous sections has been implemented using the deal.II finite element library 
described in Bangerth et al. EE). Parallelism is handled by using the MPI (Message Passing Interface) library, see 
Gabriel et al. H6l l. The subdivision and mesh distribution is done by using the p4est library from Burstedde et al. iflQll . 

The rest of this section illustrates and evaluates the performance of the above algorithm. We start in Section |4.1| 
by specifying all the numerical constants required by the algorithm. The validation of the transport code for solving 
the level-set equation in done Section [42] The validation of the two-phase fluid system is done in Section [43] 

4.1. Numerical Parameters 

Our algorithm involves several numerical parameters. In this section, we briefly recall their meaning, where 
they appear, and we specify the value of each of them. Unless specified otherwise these values are fixed for this 
entire section. A complete list of all the parameter is shown in Table [I] The parameter fi determining the width of 
the hyperbolic tangent filter ( [12] ) is defined to be fi = min^ e 7 ~ h R . This definition is used in the expression of the 
hyperbolic tangent filter ( [T2| ), in the threshold of the approximate Heaviside function ( [42] ) and the approximate sign 
function ( [38] ), and in the refinement strategy <|46]>-<|47]>. 
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Ccfl 

Ca | Ch \ Cs 

Cr | Cc | Rmax 

C\An | Ce nt | P | ^ stab 

Purpose 

CFL 

Reinitialization 

Adaptivity 

Stabilization 

Appears in 

<27)&|36| 

l 39 l 

\42\ 

l 38 l 

\46\ 

\41\ 



l 18 l 

l 17 l 

i 33 l 

Value 

0.25 

0.01 

1.25 

0.5 

2. 

2. 

2. 

0.1 

0.1 

20 

0.1 


Table 1: Description and values of the numerical parameters. 


4.2. Transport of the Lev el-Set 
4.2.1. Convergence Tests 

The consistency of the algorithm for the approximation of the level-set equation (|23])-(|24|) is a priori third-order in 
time and space. To evaluate whether this is indeed the case, we solve the linear transport equation in the unit square 
A = (0, l) 2 using the velocity field 

_/ - sin 2 (7rv) sin(27ry)cos(7zt/0.1) \ 

U t \ sin 2 (7ry) sin(27rv)cos(7r/ L /0.1) )' 

This flow is time-periodic of period 1, wich implies that 0(x, 0.2) = (po(x), for all x e A. The initial level-set function 
is chosen to be the signed distance to the line of equation y = 0.5: 


Mvy) = y- o.5. (48) 

The errors are evaluated at t = 0.2. Three different scenarios are considered: (i) no stabilization and no reinitialization; 
(ii) entropy viscosity stabilization and no reinitialization; (iii) entropy viscosity and reinitialization, i.e., the complete 
algorithm. Except for y£>, the values of the numerical parameters are given in Table [I] We consider four computations 
done on four uniform meshes with constant time steps. The mesh-size and the time step are divided by 2 each time. 
The meshes are composed of 1089, 4225, 16641and 66049 Q 2 degree of freedoms. The space discretization is chosen 
fine enough not to influence the time error. We report in Table [2j the errors \\<T> N - 0oIIl 2 (A) f° r the three scenarios and 
the observed rates of convergence. Note that the reinitialization is turned on in scenario (iii), which implies that the 
exact solution at t - 0.2 is not 0o anymore; we must instead compare with /3tanh In this particular case 

we keep ft constant to ascertain the third-order consistency of the algorithm; we set ft = 0.0203125. 



Scenario (i) 

Scenario (ii) 

Scenario (iii) 

5t 

L 2 error 

Observed rate 

L 2 error 

Observed rate 

L 2 error 

Observed rate 

le-2 

1.2e-5 

- 

1.5e-5 

- 

1.3e-3 

- 

5e-3 

1.5e-6 

3. 

1.8e-6 

3. 

1.5e-4 

3. 

2.5e-3 

1.9e-7 

3. 

2.0e-7 

3.2 

2.0e-5 

3. 

1.25e-3 

2.4e-8 

3. 

2.5e-8 

3.1 

2.6e-6 

3. 


Table 2: Single vortex test case: L 2 -norm of error and observed convergence rates for each scenario after one period 
for different time discretization resolutions, (i) RK3 Galerkin algorithm; (ii) RK3 algorithm with entropy viscosity; 
(iii) RK3 with entropy viscosity and reinitialization. Third-order convergence rate is observed in all scenarios. 


4.2.2. Rotation of a circular lev el-set: Effect of Different Viscosities and Reinitialization 

In this test case the initial data for the level set is the hyperbolic tangent filter applied to the distance to the circle 
centered at (0.5,0) and of radius 0.25 


fo(x,y) := /3 tanh 


0.25 - ((v - 0.5) 2 + y 2 )2 


(x,y)eA:=(-l,l) 2 , 


see Figure |4a| The level-set is transported by using the solid rotation velocity field: 


u (x,y) := 


-y 


(49) 
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The computation is done with the adaptive algorithm described in Section [33] using R max = 2. The initial mesh is 
uniform with hx = 0.015625 for all K e T. The resulting time-dependent mesh is such that min 7 ^ 7 - hx ~ 0.00390625, 
i.e., we set (3 = 0.00390625. The values of the other numerical constants are provided in Table [I] The time step St is 
chosen to be uniform and equal to nx 10 -3 . We compare the initial level-set with its approximation after one revolution. 
Figures [4b|[4e| illustrate the benefits of using the entropy stabilization and the reinitialization technique by comparing 
the graphs of the exact and approximate level-set functions along the line of equation x = 0.5 after one revolution. 



Figure 4: (a) Initial level-set: the dark region corresponds to 0o > 0; (b)-(e) Comparison of the graphs of the exact 
(solid line) and approximate (dashed line) level-set functions along the line of equation x = 0.5 after one revolution: 
(b) without stabilization and reinitialization (Cun = CEnt = Cx = 0 ); (c) with first-order linear stabilization, Cun = 
0.1; (d) With entropy residual stabilization, Cun = 0.1 and Ce n t = 0.1; (e) with entropy residual stabilization and 
reinitialization, Cx - 0.01. The transport without stabilization of the nearly discontinuous level-set function yields 
spurious oscillations. These oscillations are removed by the linear viscosity at the expense of large numerical diffusion. 
The numerical diffusion is minimized by turning on the entropy residual term. Finally, adding the reinitialization 
allows to nearly recover the exact profile. 


4.2.3. 3D Slotted Disk: Long Time Behavior 

A typical benchmark for the transport of a level-set function is the so-called Zalesak disk documented in IB 6 j . We 
consider in this section the three-dimensional version thereof. The computational domain is A := (0, l) 3 . The initial 
level-set is the characteristic function of a slotted sphere centered at (0.5,0.75,0.5) with a radius of 0.15. The width, 
height and depth of the slot are 0.0375, 0.15, 0.3 respectively; see Figure [5a|. The initial profile is transported by 
using the following velocity field: 


u(x, y, z) 


^ 0.5 -y ' 
x - 0.5 


v 


0 


(50) 


The time step is chosen to be uniform St = 7 rxl 0 -3 . The initial subdivision is composed of cells of diameter 0.015625. 
The adaptive mesh refinement technique described in Section 3.5 is used with R mSLX = 2; the minimum mesh size is 
minxerhx = 0.00390625. The numerical constants are given in Table [I] The computation is done until the slotted 



(a) (b) (c) (d) (e) (f) 

Figure 5: Rotating Slotted three dimensional sphere. From left to right the dark regions corresponds to O > 0 after 
0,1,2,3,4 and 5 full rotations. Oscillations and numerical diffusion are controlled by the entropy viscosity as well as 
the reinitialization algorithm. 
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sphere has undergone 5 full revolutions. The iso-surface ® = 0 is shown in Figure [5] after each of the 5 periods. 
Oscillations and numerical diffusion are controlled by the entropy viscosity and the reinitialization algorithm. A 
closer look at the slotted region is provided in Figure [6] 




Figure 6: Iso-contour ® = 0 in the plane z = 0: (A) initial data; (B) after one rotation; (C) after three rotations; (D) af¬ 
ter five rotations. Numerical diffusion is observed but is greatly minimized by the entropy viscosity and reinitialization 
techniques. 


4.2.4. Single Vortex: Large Deformations 

The Single Vortex problem consists of the deformation of a sphere by a time-periodic incompressible vortex-like 
flow. The computational domain is A = (0, l) 3 , and the time-periodic velocity field is defined by 


u (x,y 9 1) 


- sin 2 (nx) sin(27ry) cos(7z?/4) 
sin 2 (7ry) sin(27rv) cos(7zt/4) 


(51) 


The initial level-set is given by 


0o (x) := /3tanh 


dist(S,x)\ 

p r 


(52) 


where S is the sphere centered at (0.5,0.75,0.5) of radius of 0.15. The field 0o is a regularized version of the distance 
function using the tanh cut-off filter ( [T2] ). The divergence-free velocity field severely deforms the level-set until t = 2 
and returns it to its initial shape at t = 4. The time step is chosen to be uniform St = 0.001 and the final time is t = 4 
(1 cycle). The initial subdivision is made of uniform cells of diameter 0.015625; the minimum mesh size allowed 
in the adaptive mesh refinement is 0.00390625. The numerical constants are given in Table [T] Figure [7] shows the 
iso-contour ® = 0 at different times. The undeformed sphere is recovered after one cycle. 





Figure 7: 3D Vortex using From left to right: t = 0, t = l, t = 3, t = 4. The initial sphere is recovered after one 
cycle. The vectors indicate the direction of the velocity field. 


4.3. Two Phase Flows 

4.3.1. Rising Bubble: Surface Tension Benchmark 

We start the validation of the two phase flow system with the Rising bubble benchmark problem, wee e.g. Hysing 
et al. |[23l . The computational domain is A = (0, l)x(0,2), the initial data 0o is the characteristic function of a circular 
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bubble of radius 0.25 centered at (0.5,0.5). Two different sets of physical constants are considered, see Table [3j <x 
is the surface tension coefficient and g is the magnitude of the non-dimensional gravitational force. The domains 
fl + and kl~ are the domains inside and outside the bubble, respectively. The no-slip boundary condition is imposed 
at the top and bottom of the computational domain. The free slip condition ( [2c| ) is imposed on the side walls. The 


Test case 

P + 

p 


p 

g 

cr 

1 

1000 

100 

10 

1. 

0.98 

24.5 

2 

1000 

1 

10 

0.1 

0.98 

1.96 


Table 3: Two different sets of physical constants for the rising bubble benchmark problems. The first test case consists 
of a density/viscosity ratio of 10 with a large surface tension coefficient. The density/viscosity ratio is equal to 100 in 
the second case and the surface tension coefficient is smaller. 

initial subdivision is made of uniform cells of diameter 0.03125; the minimum mesh size allowed in the adaptive mesh 
refinement is 0.0078125. The time steps 6t are chosen uniform according to the CFL restriction ([27]) and the values of 
the numerical constants are given in Table[l] We compare in Figure[8]our results with those from three other methods. 
We show in panel (a) the time history of the center of mass X c := xdx/ ldx , in panel (b) the rising velocity 

u c := Uydx/ ldx , and in panel (c) the shape of the bubble at t = 3. The results of our simulations are within the 
range of those given by the benchmark algorithms. The shape of the bubble for the test case #2 at t = 3 is reported in 



(a) Center of Mass, case #1 (b) Rise Velocity, case #1 ( c ) Zero l eve l- set at t ~ 3, case #1 



Figure 8: Rising bubble test cases #1 and #2. Our simulations are within the range of the benchmark. 

Figure [9] and compared with the shapes obtained by the other algorithms described in Hysing et al. fl23l . 

4.3.2. Buckling fluids 

We now test the algorithm in the context of fluid buckling, see for instance Ville et al. M, Tome and McKee 
(321, Bonito et al. □ The test consists of letting a free-falling jet of very viscous fluid impinge on a horizontal 
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Figure 9: Rising bubble case #2 at t = 3. Different algorithms described in Hysing et al. l23l and our simulation in 
the rightmost panel. The shapes of the bubbles are qualitatively similar. The left figures are courtesy of S. Turek bub 

ffl- 


surface. The diameter of the impinging jet is 0.1 m and the inflow velocity is 1 m s -1 in A = (0,1 m) 2 . The physical 
parameters chosen for the falling fluid are p + = 1800 kg m -3 , p + = 250 Pa s and p~ = 1 kg m -3 , p~ - 2 x 10 -5 Pa s for 
the ambient fluid as in Ville et al. (35). The no-slip boundary condition (u = 0) is imposed at the bottom boundary, 
and an inflow boundary condition is imposed where |v- 0.51 < 0.05 at the top boundary. The open boundary condition 
is enforced on all the other boundaries, i.e., (2pV s u- pl)v = 0. Surface tension is neglected for this test case (cr = 0). 
The numerical constants are given in Table [T] 

The initial subdivision is composed of cells of uniform diameter 0.015 625 m, and the minimum cell diameter 
reached during the mesh adaption process is 0.003 906 m. The time evolution of the fluid is shown in Figure [lO] 
Buckling occurs after the viscous fluid impacts the rigid bottom plate. 



ii 


Figure 10: (From left to right) Time evolution of a jet of very viscous liquid inside a cavity filled with air. Buckling 
occurs when the liquid hits the floor. 



Figure 11: (a)-(e) Time evolution of a three-dimensional jet of Silicone oil falling in a cavity filled with air at time (a) 
t = 0.00501 s, (b) t = 0.0209 s, (c) t = 0.031 s, (d) t = 0.04099 s, and (e) t = 0.048998 s. For comparison, (f) shows 
the shape of the Silicone oil jet obtained in laboratory with the same physical conditions. Our numerical simulations 
are qualitatively in agreement with the physical experiment. 

We also show in Figure |TT| a three-dimensional simulation. The computational domain is A = (0,0.008 m) 3 , the 
initial diameter of the jet is 0.0004m and the inflow velocity is 1.75 ms -1 . The viscosity and density of the fluid 
are those of silicone oil: p + = 960kgm -3 , p + = 5 Pas. The viscosity and density of the ambient fluid are those of 
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air: p~ = 1.2kgm~ 3 , p~ = 2 x 10 -5 Pas. Surface tension is accounted for, cr = 0.021 Nm" 1 . The no-slip boundary 
condition (u = 0) is imposed at the bottom of the box. The inflow boundary is the disk yj(x- 0.004) 2 + (y - 0.004) 2 < 
0.0002 at the top boundary (z = 0.008). Open boundary conditions ( 2pV s u - pl)v = 0 are applied on the rest of the 
boundary. The results are obtained with adaptive mesh refinement with the minimum cell diameter 6.25 x 10~ 5 m. 
The time steps follow the CFL restriction ( [27] ) and the numerical constants are given in Table [T] 


5. Numerical Simulations of Bouncing Jets 

We now use our algorithm to predict the bouncing effect of jets of Newtonian (Section [571] ) and non-Newtonian 
(Section |52] )fluids. In both cases, the formation of a thin layer of air between the jet and the bulk of the fluid is a 
critical ingredient to observe the bouncing effect. The adaptive mesh refinement strategy adopted in our algorithm 
allows to capture this thin layer with a reasonable number of degrees of freedom. 


5.1. Two-dimensional Newtonian Bouncing Jets 

We start with a Newtonian fluid falling into a translating bath as in the experiment proposed in Thrasher et al. 
ED. The fluid is a silicone oil with viscosity p + = 0.25Pas and density p + = 960kgm 3 . The ambient fluid is air: 
p~ = 1.2kgm -3 and p~ = 2 x 10 -5 Pas. The computational domain is A := (0,0.04cm)x(0,0.06cm). The radius of 
the incoming jet is 0.25 cm; its velocity is u = (0, -5 cm s -1 ) on {(v,y) g <9A : y = 4, \x- 0.011 < 0.0025}. The region 
{(x,y) 6 A : 0 < y < 0.02} is filled by the same fluid, called the “bath”, and it moves to the right with a horizontal 
velocity Vb to be specified later. Slip boundary condition ( |2c| ) are imposed at the bottom of the cavity and the open 
boundary condition (2 pV s u - pi)v = 0 is applied to the rest of the boundary 


{(x,y) e d A : y = 4, \x - 0.011 > 0.0025} U {(x,y) e DA : x = 0 ,y> 0.02} U {(x,y) g <9A : x = 0.04}. 

The values of the numerical constants are those given in Table [T] The minimum cell diameter resulting from the 
adaptive mesh refinement strategy is 7.8125 x 10 -5 cm. The time steps are chosen to satisfy the CFL restriction ([27]). 

As already noted in Thrasher et al. ED, there is a range of velocities for which bouncing occurs, but the jet slides 
alon g the surface of the bath when the horizontal ve locity of the bath is too high. This is illustrated in Figur es [I2a| and 
|l2b[ the jet bounces when Vb = 8 cm s -1 (see Figure 12a) but it slides when V B = 25 cm s -1 (see Figure 12b). Note that 
it is necessary to include the surface tension to keep the jet stable after impinging on the free surface as illustrated in 
Figure [12c] where a simulation without surface tension is presented. Observe finally that all our numerical simulations 
show that a thin layer of air is formed between the jet and the bath each time the jet bounces. 


5.2. Kaye Effect 

The Kaye effect is the name given to the bouncing jet phenomenon when the bath is stationary. This effect has 
been observed to occur only with non-Newtonian fluids. It is now recognized that shear-thinning viscosity is a critical 
component of the Kaye effects, see e.g. Collyer and Fischer [121, Versluis et al. (34), Binder and Landig 0. Following 
Versluis et al. [34], we adopt the model of Cross m in the rest of the paper 


Mr) =Vo° + 


Po Poo 



(53) 


where po is the viscosity at zero shear stress, poo is the limiting viscosity for large stresses, y is the Frobenius norm of 
the rate-of-strain tensor 


r := 


( d 


I|V 5 U|| 


ZZ; 

Vi=1 j=\ 


duy da A 

dxi dxj) 


2 


and y c , n are two additional parameters. As a benchmark, we consider a commercial shampoo for which the shear¬ 
thinning constants corresponding to the above model ( [53] ) have been identified experimentally in Lee et al. [[25). Recall 
that viscosities cannot be measured directly but are deduced from velocity and displacement measurements so as to 
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(C) 


Figure 12: Newtonian bouncing jet (from left to right). The white arrow indicates the translation direction of the bath, 
(a) Bath velocity V B = 8 cm s -1 and surface tension coefficient cr = 21 mN m -1 ; the incoming jet bounces away from 
the bath and we observe the apparition of an air layer between the jet and the bath, (b) Bath velocity V B = 25 cm s -1 
and surface tension coefficient cr = 21 mN m _1 ; in this case, the bath velocity is too large and the jet slides along the 
bath surface; (c) Bath velocity V B = 8 cm s -1 but without surface tension cr - 0 mN m -1 ; compare with (a). 


match in some least-squares sense a behavior conjectured a-priori. The parameters yUo, p™, y c » n corresponding to the 
model ( [53] ) obtained in Lee et al. [25) are 


yuo = 5.7Pas, y Uoc = 10 3 Pas, y c - 15s \ and n = 1. (54) 

We show in Figure |T3] snapshots of experiments done with the shampoo poured at different flow rate. It was observed 
unambiguously in Lee et al. f25l that the jet slides on a lubricating air layer. 

In order to reproduce qualitatively the above experiments we consider the two-dimensional computational domain 
A = (0.496 m, 0.594 m)x(0,0.016 m). The no-slip boundary condition (u = 0) is imposed at the bottom of the 
computational domain and an inflow boundary condition is imposed on the disk {(x, y) \ \x - 0.5| < 0.00021875} on the 
top of the box. This correspond to a jet of radius 0.4375 mm. The inflow velocity is taken to be 1.75 m s -1 . The open 
boundary condition (2jjW s u - pl)v = 0 is applied on all the other boundaries. The numerical constants used in the 
simulations are listed in Table [I] The minimum mesh size attainable by adaptive mesh refinement is 3.125 x 10 -5 m, 
and time steps are chosen to comply with the CFL restriction ( [27] ). The physical parameters chosen for the fluid 
are p + = 1020 kg m -3 , the shear-thinning viscosity constants are provided in ( [55] ). We take p~ = 1.2 kg m -3 and 
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(a) 5 ml/min (b) 6 ml/min (c) 7 ml/min (d) 8 ml/min 


Figure 13: Shampoo poured at different flow rates (a) 5mLmin _1 ; (b) 6mLmin -1 ; (c) 7mLmin _1 ; and (d) 
8 mL min -1 ; No bouncing is observed at low flow rate (a) and (b). However, the Kaye effect is observed at higher flow 
rate (c) and (d). 


=2 x 10 5 Pas for the air. Surface tension is applied at the fluid/air interface with the surface tension coefficient 
cr = 0.03 N m 1 . 

After many numerical experiments it turned out that the above physical parameters did not give any Kaye ef¬ 
fect. Our interpretation is that the shear-thinning effect given by the set of parameters © is too strong; with these 
coefficients the fluid instantly becomes water-like when hitting the bottom of the cavity, see Figure |T4a| Several expla¬ 
nations for this mismatch are plausible: (i) as mentioned above the parameters ([54]) are measured under the assumption 
that the shear-thinning law follows the Cross model; (ii) the shear is over-predicted by our algorithm; (iii) the air layer 
width for this range of parameters is too thin to be captured by the algorithm; (iv) a fundamental component is missing 
in our mathematical model. These observations have lead us to consider a different set of shear-thinning parameters 
requiring a larger shear for a notable reduction of the viscosity and a smoother transition from the maximum to the 
minimum values of the viscosity, see Figure [l4b| The parameters that we now consider are the following: 

yUo = 5.7Pas, yUoo = 10 -3 Pas, y c = 970s -1 and n- 3. (55) 


This set of parameter produces the Kaye effect as shown in Figure |T5] Here again we notice that there is a very thin 



(a) Jet with parameters ( [54] ) 


Shear thinning viscosity for the Kaye Effect 



Figure 14: (Left) Falling jet of liquid described by the parameters given in ( [54] ). The shear-thinning effect is too strong 
and the fluid becomes instantly water-like when hitting the bottom of the cavity. (Right) Viscosity vs Shear for the 
parameters ( [54] ) (dotted line) and ( [55] ) (solid line). 


layer of air between the bouncing jet and the rest of the fluid, thereby adding to the large body of evidence pointing 
at the importance of air layers in bouncing jets. Let us finally mention that these computations show also that mesh 
adaptivity is critical to reproduce numerically the Kaye effect. 
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mm 


Figure 15: (Left to Right, Top to Bottom) Numerical simulation of the Kaye effect with adaptive meshes (from left to 
right and top to bottom). The viscosity parameters are given in ( [55] ). The first and last frames illustrate the adaptive 
subdivision generated by the adaptive strategy. An air layer appears in the last frame in the first row and is fully 
developed in the first frame of the second row. The apparition of the air layer coincides with the beginning of the 
bouncing phenomenon. 
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